function [u] = FEM(nelx,nely,coords,nodes,e,nu,lv,lp,phi)
%Function of FEM analysis

ne = nelx*nely; %Number of elements

% discretize the right hand side 
f = zeros(2*(nelx+1)*(nely+1), 1); %Global load vector

% if using elementStiffMatrix, get rid of *h!
f(2*lp-1:2*lp) = lv; %Substituting the local load vector into the global one

gk = grobalStiffnessMatrix(coords,ne,nodes,nelx,nely,e,nu,phi); %Construct global stiffness matrix

% distinguish different type of d.o.f.
fixedDofs = sort([2*((nelx+1)*((1:nely+1)-1)+1)-1 2*((nelx+1)*((1:nely+1)-1)+1)]); %Degenerate
allDofs = 1:2*(nelx+1)*(nely+1);
freeDofs = setdiff(allDofs,fixedDofs);

% solve the finite element equation with direct method built in matlab
u(freeDofs,:) = gk(freeDofs,freeDofs) \ f(freeDofs,:); %Solve linear system

% fixed the Dirichlet boundary d.o.f.
u(fixedDofs,:) = 0;